library(lme4)
library(stargazer)
library(performance)
library(ggplot2)
theme_set(theme_bw(base_size=12))

################################################################################
# Models using Academic age 

load("../Data/Authors.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
tsfp.hms <- lmer(model, HMS)
tsfp.ls <- lmer(model, LS)
tsfp.ps <- lmer(model, PS)
tsfp.sse <- lmer(model, SSE)

stargazer(tsfp.hms, tsfp.ls, tsfp.ps, tsfp.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0,type="text")

stargazer(tsfp.hms, tsfp.ls, tsfp.ps, tsfp.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################

# Mobility data models
load("../Data/MobilityData.RData")

mobility$week <- lubridate::week(mobility$date)
summary(mobility)
startW <- 6
stopW <- 22
cov.ag <- aggregate(cbind(workplaces,residential) ~ iso3c, data=subset(mobility, week >= startW & week <= stopW), FUN=base::mean)

HMS <- merge(HMS, cov.ag, by.x="CountryCode", by.y="iso3c", all.x=T)
LS <- merge(LS, cov.ag, by.x="CountryCode", by.y="iso3c", all.x=T)
PS <- merge(PS, cov.ag, by.x="CountryCode", by.y="iso3c", all.x=T)
SSE <- merge(SSE, cov.ag, by.x="CountryCode", by.y="iso3c", all.x=T)

model <- delta.sub ~ Women1*TimeSinceFirstPub + Women1*residential + (1 | CountryCode)
m.hms <- lmer(model, HMS)
m.ls <- lmer(model, LS)
m.ps <- lmer(model, PS)
m.sse <- lmer(model, SSE)

stargazer(m.hms, m.ls, m.ps, m.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(m.hms, m.ls, m.ps, m.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Residential","Women$\\times$Academic age", "Women$\\times$Residential","Intercept"), type="latex")


################################################################################
# Models using seniority
load("../Data/Authors.RData")

# Women1 models
nm.hms <- lmer(delta.sub ~ Women1 * Seniority + (1 | CountryCode), HMS)
nm.ls <- lmer(delta.sub ~ Women1 * Seniority + (1 | CountryCode), LS)
nm.ps <- lmer(delta.sub ~ Women1 * Seniority + (1 | CountryCode), PS)
nm.sse <- lmer(delta.sub ~ Women1 * Seniority + (1 | CountryCode), SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Doctor", "Professor", "Women$\\times$Doctor", "Women$\\times$Professor", "Intercept"), type="latex")

################################################################################
# Women2 models
model <- delta.sub ~ Women2 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0,type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################
# Research papers only models

load("../Data/AuthorsResearchPapers.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")


################################################################################
# Q1 journals only models

load("../Data/AuthorsQ1.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################
# First author only models

load("../Data/AuthorsFirst.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################
# Solo author only models

load("../Data/AuthorsSolo.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################
# COVID-only models

load("../Data/AuthorsCovid.RData")

model <- delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode)
nm.hms <- lmer(model, HMS)
nm.ls <- lmer(model, LS)
nm.ps <- lmer(model, PS)
nm.sse <- lmer(model, SSE)

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, type="text")

stargazer(nm.hms, nm.ls, nm.ps, nm.sse, report="vcsp", object.names=T, model.numbers=F, digit.separate=0, covariate.labels=c("Women", "Academic age", "Women$\\times$Academic age",  "Intercept"), type="latex")

################################################################################
# COVID-only model for HMS (Table 1)
library(lmerTest)
library(xtable)
load("../Data/AuthorsCovid.RData")

hms.only <- lmer(delta.sub ~ Women1 * TimeSinceFirstPub + (1 | CountryCode), HMS)
s <- summary(hms.only)
s

rownames(s$coefficients) <- c("Intercept", "Women", "Academic age", "Women$\\times$Academic age")
xtable(s$coefficients[,-3], digits=c(0,3,3,3,3))

s$logLik
